Fisk distribution (fisk) — the log-logistic workhorse#
The Fisk distribution (SciPy: scipy.stats.fisk) is widely known in applied statistics as the log-logistic distribution.
It is a continuous distribution on the positive real line with Pareto-like tails.
A key representation makes much of its behavior intuitive:
If (X \sim \mathrm{Fisk}(c, \text{scale}=s)) with loc=0, then
[ \log X \sim \mathrm{Logistic}(\mu=\log s,; \sigma = 1/c). ]
So you can think of Fisk as “a logistic distribution on the log-scale”.
What you’ll learn#
how to write the PDF/CDF/quantile function (and how SciPy parameterizes them)
which moments exist (and why many do not)
how to sample with NumPy only via inverse-CDF / logistic tricks
how to fit and validate the model with
scipy.stats.fisk
import platform
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
import scipy
from scipy import optimize, stats
pio.templates.default = "plotly_white"
# CKC convention for Plotly notebooks
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)
print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)
# Example parameters used throughout
c0 = 5.0
scale0 = 2.0
loc0 = 0.0
Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0
1) Title & Classification#
Name:
fisk(Fisk / log-logistic distribution; SciPy:scipy.stats.fisk)Type: Continuous
Support (standard): (x \in (0, \infty))
Parameter space (standard): shape (c>0), scale (s>0)
SciPy location/scale:
loc \in \mathbb{R},scale > 0with [ X = \mathrm{loc} + \mathrm{scale},Z,\quad Z \sim \mathrm{Fisk}(c, 1). ]
We write (2-parameter form):
[ X \sim \mathrm{Fisk}(c, s). ]
The standard form is (\mathrm{Fisk}(c,1)), i.e. stats.fisk(c, loc=0, scale=1).
2) Intuition & Motivation#
2.1 What it models#
The Fisk/log-logistic distribution models strictly positive quantities where:
values can vary over orders of magnitude, and
the right tail is heavy (polynomial decay), and
on the log-scale, the distribution looks approximately logistic.
A nice way to see the logistic link is to rewrite the CDF in terms of (y = \log x) (with loc=0):
[ F(x) = \frac{1}{1 + (s/x)^c} = \sigma\bigl(c(\log x - \log s)\bigr), \qquad \sigma(t) = \frac{1}{1+e^{-t}}. ]
So the “transition” happens around (x \approx s) (because (\log x - \log s \approx 0)).
2.2 Typical real-world use cases#
Survival analysis / reliability: time-to-event with a hazard that can rise and fall.
Economics: income/wealth-like quantities with heavy tails.
Hydrology / environmental extremes: positive heavy-tailed magnitudes.
Generative modeling: as a simple positive heavy-tailed component.
2.3 Relations to other distributions#
Log-logistic: Fisk is the log-logistic distribution.
Logistic on the log-scale: (\log X) is logistic.
Burr Type XII: Fisk is a special case of Burr XII (shape parameter (k=1)).
Pareto-like tails: (\mathbb{P}(X>x) \sim (s/x)^c) for large (x), so (c) behaves like a tail index.
3) Formal Definition#
We use the common 2-parameter form with shape (c>0) and scale (s>0).
3.1 PDF#
For (x>0):
[ f(x\mid c,s) = \frac{c}{s},\frac{(x/s)^{c-1}}{\bigl(1 + (x/s)^c\bigr)^2}. ]
3.2 CDF#
For (x>0):
[ F(x\mid c,s) = \frac{1}{1 + (s/x)^c} = \frac{(x/s)^c}{1 + (x/s)^c} = \sigma\bigl(c(\log x - \log s)\bigr). ]
The survival function is
[ \bar F(x) = 1 - F(x) = \frac{1}{1 + (x/s)^c}. ]
3.3 Quantile function (inverse CDF)#
For (u\in(0,1)):
[ Q(u) = F^{-1}(u) = s\left(\frac{u}{1-u}\right)^{1/c}. ]
3.4 SciPy parameterization#
SciPy uses:
stats.fisk(c, loc=0, scale=1)
with (z = (x-\mathrm{loc})/\mathrm{scale}) and support (z>0).
def _log1pexp(a: np.ndarray) -> np.ndarray:
"""Stable log(1+exp(a)) for real a (NumPy-only)."""
a = np.asarray(a, dtype=float)
out = np.empty_like(a, dtype=float)
pos = a > 0
out[pos] = a[pos] + np.log1p(np.exp(-a[pos]))
out[~pos] = np.log1p(np.exp(a[~pos]))
return out
def _expit(a: np.ndarray) -> np.ndarray:
"""Stable logistic sigmoid 1/(1+exp(-a)) (NumPy-only)."""
a = np.asarray(a, dtype=float)
out = np.empty_like(a, dtype=float)
pos = a >= 0
out[pos] = 1.0 / (1.0 + np.exp(-a[pos]))
ea = np.exp(a[~pos])
out[~pos] = ea / (1.0 + ea)
return out
def _logit(p: np.ndarray) -> np.ndarray:
"""Stable log(p/(1-p)) for p in (0,1)."""
p = np.asarray(p, dtype=float)
return np.log(p) - np.log1p(-p)
def fisk_logpdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Log-PDF of Fisk(c, loc, scale) in SciPy's parameterization."""
x = np.asarray(x, dtype=float)
c = float(c)
loc = float(loc)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
z = (x - loc) / scale
out = np.full_like(z, fill_value=-np.inf, dtype=float)
mask = z > 0
if not np.any(mask):
return out
logz = np.log(z[mask])
a = c * logz # log(z^c)
out[mask] = (
np.log(c)
- np.log(scale)
+ (c - 1.0) * logz
- 2.0 * _log1pexp(a)
)
return out
def fisk_pdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""PDF of Fisk(c, loc, scale)."""
return np.exp(fisk_logpdf(x, c=c, loc=loc, scale=scale))
def fisk_cdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""CDF of Fisk(c, loc, scale)."""
x = np.asarray(x, dtype=float)
c = float(c)
loc = float(loc)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
z = (x - loc) / scale
out = np.zeros_like(z, dtype=float)
mask = z > 0
if not np.any(mask):
return out
logz = np.log(z[mask])
out[mask] = _expit(c * logz)
return out
def fisk_sf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Survival function 1 - CDF of Fisk(c, loc, scale)."""
x = np.asarray(x, dtype=float)
c = float(c)
loc = float(loc)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
z = (x - loc) / scale
out = np.ones_like(z, dtype=float)
out[z <= 0] = 1.0
mask = z > 0
if not np.any(mask):
return out
logz = np.log(z[mask])
out[mask] = _expit(-c * logz)
return out
def fisk_ppf(u: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Quantile function (inverse CDF) of Fisk(c, loc, scale)."""
u = np.asarray(u, dtype=float)
c = float(c)
loc = float(loc)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
if np.any((u < 0) | (u > 1)):
raise ValueError("u must be in [0, 1]")
out = np.full_like(u, fill_value=np.nan, dtype=float)
out[u == 0] = loc
out[u == 1] = np.inf
mask = (u > 0) & (u < 1)
if not np.any(mask):
return out
out[mask] = loc + scale * np.exp(_logit(u[mask]) / c)
return out
4) Moments & Properties#
A useful fact is that the right tail is polynomial:
[ \mathbb{P}(X>x) = \bar F(x) = \frac{1}{1 + (x/s)^c} \sim (s/x)^c\quad (x\to\infty). ]
So (c) is a tail index: smaller (c) means heavier tails.
4.1 Raw moments (and existence)#
For (k>0), the raw moment exists iff (c > k), and
[ \mathbb{E}[X^k] = s^k,\Gamma!\left(1+\frac{k}{c}\right)\Gamma!\left(1-\frac{k}{c}\right) = s^k,\frac{\frac{k\pi}{c}}{\sin\left(\frac{k\pi}{c}\right)}. ]
Consequences:
Mean exists for (c>1): [\mathbb{E}[X] = s,\frac{\pi/c}{\sin(\pi/c)}.]
Variance exists for (c>2): [\mathrm{Var}(X) = s^2\left(\frac{2\pi/c}{\sin(2\pi/c)} - \left(\frac{\pi/c}{\sin(\pi/c)}\right)^2\right).]
Skewness exists for (c>3); kurtosis exists for (c>4).
4.2 Location summaries#
Median: (\mathrm{median}(X)=s) (and with SciPy
loc, the median isloc + scale).Mode (for (c>1)): [\mathrm{mode}(X)=s\left(\frac{c-1}{c+1}\right)^{1/c}.] For (c\le 1), the density is decreasing and the mode is at the lower endpoint.
4.3 MGF / characteristic function#
The MGF (M_X(t)=\mathbb{E}[e^{tX}]) does not exist for any (t>0) (polynomial tails cannot beat exponential growth).
The characteristic function (\varphi_X(t)=\mathbb{E}[e^{itX}]) exists, but it does not simplify to an elementary closed form; it’s typically computed numerically.
4.4 Entropy#
Using the logistic-on-log-scale representation, the differential entropy is
[ h(X) = 2 + \log\left(\frac{s}{c}\right)\quad \text{(nats)}. ]
(Translation by loc does not change entropy; scaling multiplies it by adding (\log(\mathrm{scale})).)
def fisk_raw_moment(k: float, c: float, scale: float = 1.0) -> float:
"""Raw moment E[X^k] for Fisk(c, scale) (loc=0).
Returns np.inf if the moment diverges (c <= k).
"""
k = float(k)
c = float(c)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
if k == 0:
return 1.0
if c <= k:
return float(np.inf)
a = np.pi * k / c
return float((scale**k) * (a / np.sin(a)))
def fisk_entropy(c: float, scale: float = 1.0) -> float:
"""Differential entropy in nats for Fisk(c, scale) (loc drops out)."""
c = float(c)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
return float(2.0 + np.log(scale / c))
def fisk_mode(c: float, loc: float = 0.0, scale: float = 1.0) -> float:
"""Mode of Fisk(c, loc, scale)."""
c = float(c)
loc = float(loc)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
if c <= 1:
return float(loc) # density decreases from the boundary
return float(loc + scale * ((c - 1.0) / (c + 1.0)) ** (1.0 / c))
def fisk_moments(c: float, loc: float = 0.0, scale: float = 1.0) -> dict:
"""Mean/variance/skewness/kurtosis (if they exist) + useful summaries."""
c = float(c)
loc = float(loc)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
m1 = fisk_raw_moment(1.0, c=c, scale=scale)
m2 = fisk_raw_moment(2.0, c=c, scale=scale)
m3 = fisk_raw_moment(3.0, c=c, scale=scale)
m4 = fisk_raw_moment(4.0, c=c, scale=scale)
mean = loc + m1 if np.isfinite(m1) else np.inf
if np.isfinite(m1) and np.isfinite(m2):
var = m2 - m1**2
else:
var = np.inf
skew = np.nan
excess_kurt = np.nan
if np.isfinite(m1) and np.isfinite(m2) and np.isfinite(m3) and var > 0:
mu3 = m3 - 3 * m2 * m1 + 2 * (m1**3)
skew = mu3 / (var ** 1.5)
if np.isfinite(m1) and np.isfinite(m2) and np.isfinite(m3) and np.isfinite(m4) and var > 0:
mu4 = m4 - 4 * m3 * m1 + 6 * m2 * (m1**2) - 3 * (m1**4)
excess_kurt = mu4 / (var**2) - 3.0
median = loc + scale
return {
"mean": float(mean),
"var": float(var),
"skew": float(skew) if np.isfinite(skew) else np.nan,
"excess_kurtosis": float(excess_kurt) if np.isfinite(excess_kurt) else np.nan,
"median": float(median),
"mode": fisk_mode(c=c, loc=loc, scale=scale),
"entropy": fisk_entropy(c=c, scale=scale),
}
m0 = fisk_moments(c0, loc=loc0, scale=scale0)
m0
{'mean': 2.13791866423119,
'var': 0.7145293838425228,
'skew': 2.4852755496867136,
'excess_kurtosis': 26.556191909249108,
'median': 2.0,
'mode': 1.8442158229634555,
'entropy': 1.083709268125845}
5) Parameter Interpretation#
5.1 Scale s#
scale = ssets the median: (\mathrm{median}(X)=s) (whenloc=0).On the log-scale, it is a location parameter: (\mu = \log s).
5.2 Shape c#
ccontrols the steepness around the median and the tail heaviness.The tail index is
c: (\mathbb{P}(X>x) \sim (s/x)^c).Moments exist only up to order
< c.On the log-scale,
cis the inverse of the logistic scale: (\sigma = 1/c).
5.3 Hazard shape (survival analysis intuition)#
The hazard rate (h(x)=f(x)/(1-F(x))) (for loc=0) simplifies to
[ h(x) = \frac{c}{s},\frac{(x/s)^{c-1}}{1 + (x/s)^c}. ]
For (c \le 1), the hazard is decreasing.
For (c > 1), the hazard is typically unimodal (rises, then falls).
def fisk_hazard(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Hazard h(x) = f(x) / (1 - F(x))."""
f = fisk_pdf(x, c=c, loc=loc, scale=scale)
s = fisk_sf(x, c=c, loc=loc, scale=scale)
out = np.full_like(f, fill_value=np.nan, dtype=float)
mask = s > 0
out[mask] = f[mask] / s[mask]
return out
x = np.logspace(-3, 2, 800) * scale0
c_values = [0.7, 1.0, 2.0, 5.0]
fig = make_subplots(
rows=1,
cols=3,
subplot_titles=["PDF (log x)", "CDF (log x)", "Hazard (log x)"],
)
for c in c_values:
fig.add_trace(
go.Scatter(x=x, y=fisk_pdf(x, c=c, scale=scale0), name=f"c={c}", mode="lines"),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=x,
y=fisk_cdf(x, c=c, scale=scale0),
name=f"c={c}",
mode="lines",
showlegend=False,
),
row=1,
col=2,
)
fig.add_trace(
go.Scatter(
x=x,
y=fisk_hazard(x, c=c, scale=scale0),
name=f"c={c}",
mode="lines",
showlegend=False,
),
row=1,
col=3,
)
for j in [1, 2, 3]:
fig.update_xaxes(type="log", title_text="x", row=1, col=j)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_yaxes(title_text="h(x)", row=1, col=3)
fig.update_layout(title="Fisk distribution shape effects (scale fixed)", height=420)
fig.show()
6) Derivations#
6.1 Expectation (raw moments)#
Start from the definition for (k>0):
[ \mathbb{E}[X^k] = \int_0^\infty x^k,\frac{c}{s},\frac{(x/s)^{c-1}}{(1 + (x/s)^c)^2},dx. ]
Use the substitution (t = (x/s)^c), so (x = s,t^{1/c}) and (dx = \frac{s}{c} t^{1/c-1},dt). A pleasant simplification happens:
[ \frac{c}{s}(x/s)^{c-1},dx = \frac{c}{s},t^{(c-1)/c},\frac{s}{c}t^{1/c-1}dt = dt. ]
So
[ \mathbb{E}[X^k] = s^k \int_0^\infty \frac{t^{k/c}}{(1+t)^2},dt. ]
Recognize the Beta-function integral:
[ \int_0^\infty \frac{t^{p-1}}{(1+t)^{m}},dt = B(p, m-p),\quad 0<p<m. ]
Here (m=2) and (p=1+k/c), so the condition becomes (c>k). Then
[ \mathbb{E}[X^k] = s^k,B\left(1+\frac{k}{c},,1-\frac{k}{c}\right) = s^k,\Gamma\left(1+\frac{k}{c}\right)\Gamma\left(1-\frac{k}{c}\right). ]
Using (\Gamma(1+z)\Gamma(1-z)=\frac{\pi z}{\sin(\pi z)}) gives the sine form.
6.2 Variance#
When (c>2),
[ \mathrm{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2, ]
with (\mathbb{E}[X]) requiring (c>1).
6.3 Likelihood (i.i.d. sample)#
For i.i.d. observations (x_1,\dots,x_n > 0) (assuming loc=0), the log-likelihood is
[ \ell(c,s) = \sum_{i=1}^n \log f(x_i\mid c,s) = n\log c - n\log s + (c-1)\sum_i\log\left(\frac{x_i}{s}\right) - 2\sum_i \log\left(1+\left(\frac{x_i}{s}\right)^c\right). ]
There is no closed-form MLE for ((c,s)); we maximize (\ell) numerically. A common trick is to optimize over ((\log c, \log s)) to enforce positivity.
def fisk_loglikelihood(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> float:
"""Total log-likelihood for i.i.d. data under Fisk(c, loc, scale)."""
return float(np.sum(fisk_logpdf(x, c=c, loc=loc, scale=scale)))
def fisk_mle_loc_fixed(x: np.ndarray, loc: float = 0.0) -> tuple[float, float]:
"""Numerical MLE for (c, scale) with loc fixed.
Uses a logistic-on-log-scale initialization.
"""
x = np.asarray(x, dtype=float)
if np.any(x <= loc):
raise ValueError("all observations must be > loc")
x_pos = x - loc
# Initialization from log-scale logistic approximation
y = np.log(x_pos)
mu_init = float(np.median(y))
q25, q75 = np.quantile(y, [0.25, 0.75])
iqr = float(q75 - q25)
# Logistic IQR = 2*s*log(3), with s = 1/c
c_init = float(max(2.0 * np.log(3.0) / max(iqr, 1e-6), 1e-3))
scale_init = float(max(np.exp(mu_init), 1e-6))
def nll(theta: np.ndarray) -> float:
log_c, log_scale = float(theta[0]), float(theta[1])
c = float(np.exp(log_c))
scale = float(np.exp(log_scale))
return -fisk_loglikelihood(x, c=c, loc=loc, scale=scale)
res = optimize.minimize(
nll,
x0=np.array([np.log(c_init), np.log(scale_init)]),
method="Nelder-Mead",
options={"maxiter": 5000},
)
log_c_hat, log_scale_hat = res.x
return float(np.exp(log_c_hat)), float(np.exp(log_scale_hat))
# Quick demo on synthetic data
x_demo = stats.fisk.rvs(c0, loc=loc0, scale=scale0, size=2000, random_state=rng)
c_hat, scale_hat = fisk_mle_loc_fixed(x_demo, loc=loc0)
print("true (c, scale) =", (c0, scale0))
print("MLE (c, scale) =", (c_hat, scale_hat))
true (c, scale) = (5.0, 2.0)
MLE (c, scale) = (5.081944016481647, 1.9992827793076853)
7) Sampling & Simulation#
Inverse-transform sampling uses the quantile function.
If (U\sim\mathrm{Uniform}(0,1)), then
[ X = Q(U) = \mathrm{loc} + \mathrm{scale},\exp\left(\frac{\mathrm{logit}(U)}{c}\right) = \mathrm{loc} + \mathrm{scale}\left(\frac{U}{1-U}\right)^{1/c}. ]
Algorithm (NumPy-only):
Draw (u_1,\dots,u_n) i.i.d. from
Uniform(0,1).Clip away from exactly 0 or 1 (to avoid
±infafterlogit).Compute
z = exp(logit(u) / c)and returnloc + scale * z.
Because the distribution is logistic on the log-scale, this is equivalent to sampling
[ Y \sim \mathrm{Logistic}(\mu=\log(\mathrm{scale}), \sigma = 1/c),\qquad X = \mathrm{loc}+e^Y. ]
def fisk_rvs_numpy(
c: float,
size: int,
loc: float = 0.0,
scale: float = 1.0,
rng: np.random.Generator | None = None,
) -> np.ndarray:
"""Sample Fisk(c, loc, scale) using NumPy only (inverse CDF)."""
c = float(c)
loc = float(loc)
scale = float(scale)
if c <= 0:
raise ValueError("c must be > 0")
if scale <= 0:
raise ValueError("scale must be > 0")
if rng is None:
rng = np.random.default_rng()
u = rng.random(size)
# Avoid u=0 or u=1 exactly (would map to +/-inf)
eps = np.finfo(float).eps
u = np.clip(u, eps, 1.0 - eps)
z = np.exp(_logit(u) / c)
return loc + scale * z
# Monte Carlo validation (choose c>2 so mean/variance exist)
n = 200_000
samples_numpy = fisk_rvs_numpy(c0, size=n, loc=loc0, scale=scale0, rng=rng)
mc_mean = samples_numpy.mean()
mc_var = samples_numpy.var(ddof=0)
print("theory mean", m0["mean"], "MC", mc_mean)
print("theory var ", m0["var"], "MC", mc_var)
# Compare NumPy-only sampler to SciPy sampler (quick 2-sample KS test)
samples_scipy = stats.fisk.rvs(c0, loc=loc0, scale=scale0, size=n, random_state=rng)
ks = stats.ks_2samp(samples_numpy[:30_000], samples_scipy[:30_000])
ks
theory mean 2.13791866423119 MC 2.137557641110836
theory var 0.7145293838425228 MC 0.7035045300612408
KstestResult(statistic=0.009700000000000042, pvalue=0.11809773491898767, statistic_location=2.1769456563992517, statistic_sign=-1)
8) Visualization#
We’ll visualize:
the PDF with a Monte Carlo histogram overlay
the CDF vs an empirical CDF
the log-scale view ((\log X) should look logistic)
Because Fisk can be heavy-tailed, a log-x axis is often the clearest.
# PDF + histogram + CDF/ECDF
x_grid = np.logspace(-3, 2, 600) * scale0
# Empirical CDF (subsample for plotting)
sub = np.sort(samples_numpy[:10_000])
ecdf_y = np.arange(1, sub.size + 1) / sub.size
fig = make_subplots(
rows=1,
cols=2,
subplot_titles=["PDF with Monte Carlo histogram", "CDF vs empirical CDF"],
)
# Histogram (density)
fig.add_trace(
go.Histogram(
x=samples_numpy,
histnorm="probability density",
nbinsx=120,
name="samples",
opacity=0.5,
),
row=1,
col=1,
)
# Theoretical PDF
fig.add_trace(
go.Scatter(
x=x_grid,
y=fisk_pdf(x_grid, c=c0, loc=loc0, scale=scale0),
mode="lines",
name="theory pdf",
line=dict(width=2),
),
row=1,
col=1,
)
# Theoretical CDF
fig.add_trace(
go.Scatter(
x=x_grid,
y=fisk_cdf(x_grid, c=c0, loc=loc0, scale=scale0),
mode="lines",
name="theory cdf",
line=dict(width=2),
showlegend=False,
),
row=1,
col=2,
)
# Empirical CDF
fig.add_trace(
go.Scatter(
x=sub,
y=ecdf_y,
mode="lines",
name="empirical cdf",
showlegend=False,
),
row=1,
col=2,
)
fig.update_xaxes(type="log", title_text="x", row=1, col=1)
fig.update_xaxes(type="log", title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(height=420, title=f"Fisk(c={c0}, scale={scale0})")
fig.show()
# Log-scale view: log(X) should look logistic
y = np.log(samples_numpy - loc0)
mu_hat = np.median(y)
# logistic scale estimate via IQR: IQR = 2*s*log(3)
q25, q75 = np.quantile(y, [0.25, 0.75])
s_hat = (q75 - q25) / (2 * np.log(3))
# Compare to SciPy logistic on log-scale
logistic_dist = stats.logistic(loc=mu_hat, scale=s_hat)
y_grid = np.linspace(np.quantile(y, 0.01), np.quantile(y, 0.99), 500)
fig2 = go.Figure()
fig2.add_trace(
go.Histogram(x=y, histnorm="probability density", nbinsx=120, name="log samples", opacity=0.6)
)
fig2.add_trace(go.Scatter(x=y_grid, y=logistic_dist.pdf(y_grid), mode="lines", name="logistic fit"))
fig2.update_layout(
title="On the log-scale, Fisk becomes logistic",
xaxis_title="y = log(x)",
yaxis_title="density",
height=380,
)
fig2.show()
9) SciPy Integration (scipy.stats.fisk)#
SciPy parameterization:
stats.fisk(c, loc=0, scale=1)
cis the shape.locshifts the support to(loc, ∞).scalerescales (and in the 2-parameter form withloc=0, it equals the median).
Key methods:
pdf(x),logpdf(x),cdf(x),ppf(q)rvs(size, random_state=...)fit(data, ...)for MLE
dist = stats.fisk(c0, loc=loc0, scale=scale0)
x_test = np.array([0.5, 1.0, 2.0, 4.0, 8.0])
pdf_scipy = dist.pdf(x_test)
cdf_scipy = dist.cdf(x_test)
pdf_ours = fisk_pdf(x_test, c=c0, loc=loc0, scale=scale0)
cdf_ours = fisk_cdf(x_test, c=c0, loc=loc0, scale=scale0)
print("pdf close?", np.allclose(pdf_scipy, pdf_ours))
print("cdf close?", np.allclose(cdf_scipy, cdf_ours))
# Fitting (MLE) with SciPy; fix loc if you know the support starts at 0.
c_fit, loc_fit, scale_fit = stats.fisk.fit(samples_numpy[:30_000], floc=0)
print("fit (c, loc, scale) =", (c_fit, loc_fit, scale_fit))
pdf close? True
cdf close? True
fit (c, loc, scale) = (4.998041440572122, 0, 2.007897677685366)
10) Statistical Use Cases#
10.1 Hypothesis testing / goodness-of-fit#
Typical workflow:
Fit the distribution parameters (MLE).
Check fit with diagnostics:
KS test (quick but less sensitive in the tails)
QQ/PP plots (especially important for tail fit)
tail-specific checks (log-log plots, exceedance probabilities)
You can also compare candidates (e.g., Fisk vs lognormal vs Weibull) using AIC.
10.2 Bayesian modeling#
There is no conjugate prior for ((c,s)), but Bayesian inference is straightforward with MCMC/VI. A useful modeling trick is to work with (y_i = \log x_i):
[ y_i \sim \mathrm{Logistic}(\mu=\log s,,\sigma=1/c), ]
and add the Jacobian term if you write the model on (x) directly.
Below we’ll do a simple grid posterior approximation with priors on (\log c) and (\log s).
10.3 Generative modeling#
Because sampling is cheap, Fisk is convenient for generating positive heavy-tailed synthetic data. For example, in an accelerated-failure-time style model, you can let the scale depend on covariates:
[ \log X = \beta^\top x + \varepsilon,\quad \varepsilon \sim \mathrm{Logistic}(0, 1/c). ]
# Hypothesis testing + model comparison (AIC) on synthetic data
true_c, true_scale = 2.5, 1.7
x_data = fisk_rvs_numpy(true_c, size=1500, scale=true_scale, rng=rng)
# Fit Fisk
c_hat, loc_hat, scale_hat = stats.fisk.fit(x_data, floc=0)
ll_fisk = float(np.sum(stats.fisk.logpdf(x_data, c_hat, loc=loc_hat, scale=scale_hat)))
aic_fisk = 2 * 2 - 2 * ll_fisk # 2 params (c, scale) since loc fixed
# Fit lognormal
s_logn, loc_logn, scale_logn = stats.lognorm.fit(x_data, floc=0)
ll_logn = float(np.sum(stats.lognorm.logpdf(x_data, s_logn, loc=loc_logn, scale=scale_logn)))
aic_logn = 2 * 2 - 2 * ll_logn
# Fit Weibull
c_w, loc_w, scale_w = stats.weibull_min.fit(x_data, floc=0)
ll_w = float(np.sum(stats.weibull_min.logpdf(x_data, c_w, loc=loc_w, scale=scale_w)))
aic_w = 2 * 2 - 2 * ll_w
print("Fitted Fisk (c, scale)", (c_hat, scale_hat), "AIC", aic_fisk)
print("Fitted lognorm(s, scale)", (s_logn, scale_logn), "AIC", aic_logn)
print("Fitted Weibull(c, scale)", (c_w, scale_w), "AIC", aic_w)
# KS test against fitted Fisk CDF (note: p-values are approximate after fitting)
dist_hat = stats.fisk(c_hat, loc=0, scale=scale_hat)
ks = stats.kstest(x_data, dist_hat.cdf)
ks
Fitted Fisk (c, scale) (2.459344430482952, 1.6658920422291383) AIC 4812.347844437876
Fitted lognorm(s, scale) (0.7263480358234555, 1.661320104043914) AIC 4824.475216282033
Fitted Weibull(c, scale) (1.3345162113961129, 2.3881880476861284) AIC 5114.144353642654
KstestResult(statistic=0.01612416741026612, pvalue=0.8241228923230904, statistic_location=2.6508596780817344, statistic_sign=-1)
# Simple Bayesian grid posterior over (c, scale) with log-normal-ish priors
x_small = x_data[:200]
logx = np.log(x_small)
c_grid = np.linspace(0.8, 6.0, 120)
scale_grid = np.linspace(0.3, 4.0, 140)
log_scale_grid = np.log(scale_grid)
# Priors: log c ~ N(log 2, 0.7^2), log scale ~ N(log 2, 0.7^2)
mu0 = np.log(2.0)
s0 = 0.7
log_prior_c = stats.norm.logpdf(np.log(c_grid), loc=mu0, scale=s0)
log_prior_scale = stats.norm.logpdf(log_scale_grid, loc=mu0, scale=s0)
log_post = np.empty((c_grid.size, scale_grid.size), dtype=float)
for i, c in enumerate(c_grid):
# Vectorized over scale_grid
logz = logx[:, None] - log_scale_grid[None, :]
a = c * logz
# logpdf for each observation and scale
logpdf = (
np.log(c)
- log_scale_grid[None, :]
+ (c - 1.0) * logz
- 2.0 * _log1pexp(a)
)
loglik = logpdf.sum(axis=0)
log_post[i, :] = loglik + log_prior_c[i] + log_prior_scale
# Normalize to probabilities
log_post -= np.max(log_post)
post = np.exp(log_post)
post /= post.sum()
# Posterior summaries
c_mean = float((post.sum(axis=1) * c_grid).sum())
scale_mean = float((post.sum(axis=0) * scale_grid).sum())
ij_map = np.unravel_index(np.argmax(post), post.shape)
c_map = float(c_grid[ij_map[0]])
scale_map = float(scale_grid[ij_map[1]])
print("posterior mean (c, scale) =", (c_mean, scale_mean))
print("MAP (c, scale) =", (c_map, scale_map))
fig = go.Figure(
data=go.Heatmap(
x=scale_grid,
y=c_grid,
z=post,
colorscale="Viridis",
colorbar=dict(title="posterior"),
)
)
fig.update_layout(
title="Grid posterior p(c, scale | data) (loc fixed to 0)",
xaxis_title="scale",
yaxis_title="c",
height=420,
)
fig.show()
posterior mean (c, scale) = (2.5669389607324895, 1.697934398969099)
MAP (c, scale) = (2.547899159663866, 1.6841726618705037)
# Generative modeling sketch: scale depends on covariates (AFT-style)
n_gen = 1500
x_feat = rng.normal(size=n_gen)
beta0, beta1 = 0.2, 0.9
scale_x = np.exp(beta0 + beta1 * x_feat) # positive, log-linear
u = rng.random(size=n_gen)
u = np.clip(u, np.finfo(float).eps, 1.0 - np.finfo(float).eps)
# Same shape c0, but varying scale per observation
x_gen = scale_x * np.exp(_logit(u) / c0)
fig = go.Figure()
fig.add_trace(
go.Scatter(x=x_feat, y=x_gen, mode="markers", marker=dict(size=4, opacity=0.5))
)
fig.update_layout(
title="Synthetic data: Fisk noise with covariate-dependent scale",
xaxis_title="feature x",
yaxis_title="generated positive outcome",
height=380,
)
fig.show()
11) Pitfalls#
Invalid parameters:
c <= 0orscale <= 0is not valid.Moment non-existence:
if
c <= 1, the mean is infinite/undefinedif
c <= 2, the variance is infinitemore generally, (\mathbb{E}[X^k]) exists only for
k < c
Parameterization confusion:
“Fisk” and “log-logistic” are the same distribution.
Many texts use (\alpha) (shape) and (\beta) (scale); SciPy uses
candscale.On the log-scale, the logistic scale is
1/c.
Numerical overflow:
direct computation of ((x/s)^c) can overflow for large
x.prefer
logpdf, and compute the CDF viasigmoid(c * log(x/s)).
Sampling edge cases:
inverse CDF uses
logit(u); clipuaway from 0 and 1.
Fitting with
locfree:estimating
locchanges the support boundary and can be unstable; fixlocwhen it is known (e.g., 0 for strictly positive data).
# Numerical stability demo: naive power vs log-space
def fisk_pdf_naive(x: np.ndarray, c: float, scale: float = 1.0) -> np.ndarray:
x = np.asarray(x, dtype=float)
z = x / scale
return (c / scale) * (z ** (c - 1.0)) / (1.0 + z**c) ** 2
x_big = np.array([1e2, 1e10, 1e50])
c_test = 5.0
with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
naive = fisk_pdf_naive(x_big, c=c_test, scale=1.0)
stable = fisk_pdf(x_big, c=c_test, scale=1.0)
log_stable = fisk_logpdf(x_big, c=c_test, scale=1.0)
print("x", x_big)
print("naive pdf", naive)
print("stable pdf", stable)
print("logpdf ", log_stable)
x [1.e+02 1.e+10 1.e+50]
naive pdf [0. 0. 0.]
stable pdf [0. 0. 0.]
logpdf [ -26.0216 -136.5457 -689.1661]
12) Summary#
fiskis the log-logistic distribution: (\log X) is logistic.It models positive heavy-tailed data; the tail index is
c.Moments exist only up to order
< c(mean requiresc>1, variance requiresc>2).Sampling is easy via the inverse CDF:
x = loc + scale * exp(logit(u)/c).SciPy provides a robust implementation:
scipy.stats.fisk(pdf/cdf/rvs/fit).
References
SciPy:
scipy.stats.fiskBurr Type XII distribution (Fisk as a special case)
Logistic distribution entropy and quantiles (for the log-scale view)